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We present here our study of the adiabatic quantum dynamics of a random Ising chain across its 
quantum critical point. The model investigated is an Ising chain in a transverse field with disorder 
present both in the exchange coupling and in the transverse field. The transverse field term is 
proportional to a function Y(t) which, as in the Kibble-Zurek mechanism, is linearly reduced to 
zero in time with a rate t , T(t) = —t/r, starting at t = — oo from the quantum disordered phase 
(r = oo) and ending at t — in the classical ferromagnetic phase (T = 0). We first analyze the 
distribution of the gaps, occurring at the critical point F c = 1, which are relevant for breaking 
the adiabaticity of the dynamics. We then present extensive numerical simulations for the residual 
energy E Tas and density of defects pk at the end of the annealing, as a function of the annealing 
inverse rate r. Both the average E les (r) and pk(r) are found to behave logarithmically for large r, 
but with different exponents, [E TCB (r) / L] av ~ l/ln^(r) with £ w 3.4, and [pfc(r)] av ~ l/ln 2 (r). We 
propose a mechanism for 1/ In 2 r-behavior of [pk]nv based on the Landau-Zener tunneling theory 
and on a Fisher's type real-space renormalization group analysis of the relevant gaps. The model 
proposed shows therefore a paradigmatic example of how an adiabatic quantum computation can 
become very slow when disorder is at play, even in absence of any source of frustration. 
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I. INTRODUCTION 

How effective is it to execute a given computational 
task by slowly varying in time the Hamiltonian of a quan- 
tum system? Is it possible to find the ground state of a 
classical system by slowly annealing away its quantum 
fluctuations? What is the density of defects left over 
after a passage through a continuous (quantum) phase 
transition? These seemingly different problems, rubri- 
cated under the names of Adiabatic quantum compu- 
tation 1 -, Quantum Annealing 2 -" 3 -^^ and Kibble-Zurek 
topological defect formation in quantum phase transi- 
tions^, are ultimately related to the understanding of the 
failure of the adiabatic approximation in a many-body 
system passing through a quantum critical point. 

Adiabatic quantum computation (AQC), alias Quan- 
tum Annealing (QA), is a possible alternative to the stan- 
dard circuit-theory approach to Quantum Computation 
(QC) 9 . Indeed, as shown by Aharonov et a/. 10 ' 11 , any 
quantum algorithm can be equivalently reformulated in 
terms of the adiabatic evolution of an appropriate time- 
dependent Hamiltonian H(t) = [I — f(t)]H- m + f(t)Hf\ n , 
f(t) being a generic function of time such that /(0) = 
and /(tfin) = 1- The initial Hamiltonian H lni for which 
we know the ground state, provides the input of the algo- 
rithm. The final Hamiltonian Hj\ n is constructed appro- 
priately so as to possess the solution of the computational 



task as its ground state. The knowledge of the equiva- 
lence of computational power between the two different 
QC schemes, however, does not provide a practical way 
of constructing H m and Ha n for a given computational 
problem. Understanding what computational problems 
can be efficiently solved by AQC-QA is, in general, a very 
difficult problem. In order to solve the task one has to 
find a suitable path in Hamiltonian space in such a way 
that the resulting Schrodinger evolution efficiently drives 
the system from some simple initial quantum state |^i n ) 
to the sought final ground stat o 12 ' 13 . The accuracy of the 
computation, which relies on the possibility for the sys- 
tem to remain in the instantaneous ground state during 
the dynamics, is ultimately limited by the fact that at 
specific times the instantaneous Hamiltonian presents a 
gap between the ground and the first excited state which 
closes on increasing the size of the input. 

On totally independent grounds, the study of topolog- 
ical defect formation goes back to the 80's, motivated by 
the effort to understand signatures of phase transitions 
which have occurred in the early univers o 14 ' 15 by deter- 
mining the density of defects left in the broken symmetry 
phase as a function of the rate of quench. By means of the 
so called Kibble-Zurek mechanism, a scaling law relates 
the density of defects to the rate of quench. The sugges- 
tion of Zurek to simulate transitions in the early universe 
by means of condensed matter system has stimulated an 
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intense experimental activit y 16 ' 17 aimed at verifying the 
Kibble-Zurek theory. The understanding of defect for- 
mation was later explored also in the case of a quantum 
phase transitions^, where the crossing of the critical 
point is done by varying a parameter in the Hamilto- 
nian. These works have stimulated an intense activity 
where several different quantum systems undergoing a 
quantum phase transition were scrutinized. In the past 
couple of years there have been a number of results ob- 
tained in the area of adiabatic dynamics of many-body 
system s 19 ' 20 . Most of the works concentrated on the one- 
dimensional Ising model. Soon after the appearance of 
Ref. |8|, Dziarmag a 2 ! obtained analytically the scaling law 
for the density of defects by resorting to the exact solu- 
tion of Pfeuty 2 ^. A detailed analysis a ' la Landau-Zener 
was presented in Refs.— i 24 i 25 . The effect of an exter- 
nal noise on the adiabatic evolution and its consequences 
for the Kibble-Zurek mechanism has been discussed in2£. 
Recently, quenches in Bosc-Hubbard models were ana- 
lyze d 27 ' 28 as well. Observables which were analyzed to 
quantify the loss of adiabaticity in the critical region were 
typically the density of defects left behind in the broken 
symmetry phase, the fidelity of the evolved state with 
respect to the ground state, and, in few cases, also the 
residual block entrop y 29 ' 30 . This brief overview of recent 
works accounts only for papers dealing with adiabatic 
dynamics, without touching the vast literature treating 
the case of sudden quenches. 

In the present work we analyze the adiabatic dynamics 
in a one-dimensional quantum disordered Ising model in 
a random transverse field. The reasons for considering 
this problem are various. First of all it is an important 
ground test for the Kibble-Zurek mechanism. In addi- 
tion, although in a very simplified manner, it may help 
in understanding more interesting problems that can be 
formulated in terms of interacting Ising spins, Traveling 
Salesman 3 -^ and Satisfiability 32 problems being only two 
well-known examples. The simplicity of our test prob- 
lem lies in the particularly simple geometry of the inter- 
actions, which forbids frustration. The only ingredient 
that our problem shares with more challenging computa- 
tional tasks is the fact that the interactions are chosen to 
be random. This feature, the presence of disorder, makes 
the problem interesting and non-trivial for a physically 
inspired computational approach based on AQC-QA. 

Of particular relevance for us is RefHH where this 
model was analyzed first, and the anomalously slow 
dynamics characterized by an average density of kinks 
which vanishes only logarithmically with the annealing 
rate. Here we extend this work by presenting a detailed 
analysis of the statistics of both the residual energy and 
kink density. In a disordered chain, the formation of the 
kinks is no longer translational invariant and therefore it 
affects in a non-trivial way, as we will show below, the 
scaling of the residual energy. 

The rest of the paper is organized as follows: In Sec.HTl 
we define the problem and the technique to solve the adi- 
abatic dynamics of the random Ising chain, and next, in 



Sec. IIII1 we introduce the quantities — residual energy 
and density of defects — that we calculate to quantify 
the departure from the adiabatic ground state. In Sec.HVl 
we present our numerical results for both these quanti- 
ties, together with an analysis of the largc-anncaling-time 
behavior of the density of defects, based on the Landau- 
Zener theory, explicitly showing the slow dynamics which 
the disorder entails. In the final section we present a crit- 
ical assessment of our findings, and a concluding discus- 
sion. 



II. THE MODEL 

As discussed in the Introduction, our aim is to ana- 
lyze the adiabatic dynamics of a one-dimensional random 
Ising model defined by the Hamiltonian 

H(t) = - Ji^t+i - r(t) £ hiaf . (i) 

i i 

In the previous expression o~f {a = x, z) are Pauli matri- 
ces for the i-th spin of the chain, Jj are random couplings 
between neighboring spins, and hi are random transverse 
fields. The time-dependent function T(t) rescaling the 
transverse field term allows us to drive the system form 
a region of infinitely high transverse fields (r = oo, where 
the ground state has all spins aligned along x, see below), 
to the case of a classical Ising model (r = 0). Specifically, 
we will take in the following T(t) to be a linear function 
of time characterized by an annealing rate r _1 

T(<) = for t G (-oo,0] . 

r 

In one-dimension, and for nearest- neighbor couplings, 
there is no frustration associated to the random nature 
of the couplings Jf. by appropriately performing spin ro- 
tations of 7r along the x-spin axis, we can always change 
the desired a\ into — <r| and invert accordingly the signs 
of the couplings in such a way that all Ji's turn out to be 
non-negative. We therefore assume that the Ji are ran- 
domly distributed in the interval [0, 1], specifically with 
a flat distribution ir[J] = 0(J)6(\ — J), where 9 is the 
Heaviside function. The same distribution is used for the 
random field ir[h] = 8(h)9(l — h). This is different from 
the model considered in RefJH, where the disorder was 
introduced in the exchange coupling only. We find the 
present choice quite convenient since, by duality argu- 
ments 3 ^, the critical point separating the large-r quan- 
tum paramagnetic phase from the low-r ferromagnetic 
region is known to be located at T c = 1. 

At the initial time t- m — — oo the ground state of 
H(t- m ), completely dominated by the transverse field 
term, is simply the state with all spins aligned along the 
+x spin direction: |# in ) = J], \x)i = Ui[\1)t + \i)i\/V2. 
On the other side of the transition point T c , the final 
Hamiltonian H(t& n ) — H c i describes a random ferromag- 
net whose ground states, which we aim to reach by adia- 
batically switching off T(t), are obviously the two trivial 
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states |* T ) = njT)i and = ILU>< : as an op- 

timization problem, H^ n represents, therefore, a trivial 
problem. 

Even if the ground states in the two limiting cases, 
r = oo and Y = 0, are very easy to find, when it 
comes to dynamics, the evolution dictated by H(t) is 
no longer a trivial problem. The instantaneous spectrum 
of the Hamiltonian H(t) is gapless in the thermodynamic 
limit^. This implies that, during the adiabatic evolution, 
defects in the form of domain walls between differently 
aligned ferromagnetic ground states, of the type 

i-.-uuummmiu---) 

are formed, and reflected in a whole structure of closing 
gaps will appear in the instantaneous spectrum. 



A. Fermion representation and Bogoliubov-de 
Gennes equations 

By means of the Jordan- Wigner transformation, the 
one-dimensional Ising model is reduced to a free fermion 
model. One first writes the spin operators in terms of 
hard-core bosons and at in a representation that maps 
the state \a\ = +1) -> |1) 4 = oj|0)j and \a\ = -1) -> 
|0)i, with the hard-core constraint (a|) 2 |0)i = 0: of = 
2a\a,i — 1, af = a« + a- , and of = —i(a\ — The hard- 
core boson operators et^ are then re-expressed in terms of 

spinless fermions operators c,: a« = e l7r ^'<* c ' Ci Ci. After 
a 7r/2 rotation around the y-axis, which maps a x — > a z 
and cr 2 — * —(J x , the Hamiltonian in Eq.([T]) can be rewrit- 
ten in terms of fermion operators as 

L-l L 

H=-J2 JiicUi+i+ch+i+K-cl-W^hich > ( 2 ) 

i i 

where we have assumed open boundary conditions (OBC) 
for the spin-chain. For the case of periodic boundary 
conditions (PBC) on the spins, <jl+i = cij extra bound- 
ary terms appear in the fermionic Hamiltonian, of the 
form A#pbc = Jl{-^) Nf {c\c\ + c\c x + H.c.}, where 

Np = J^i c l c i is t ne total number of fermions. Notice 
that although Np is not conserved by the Hamiltonian 
©, the parity of Np is conserved: (—1) Nf is a constant 
of motion with value 1 or — 1. 



1. Statics 

The model in Eq. @ can be diagonalized through a Bo- 
goliubov rotatio n 35 ' 36 , by introducing the new fermionic 



operators 7 M and 7? 

L 
L 

a = ^2(ui^ + v^l) , (3) 

where the L-dimensional vectors u M and v M , for fi = 
1, • • • ,L, satisfy the Bogoliubov-de Gennes equations: 

A ■ u M + B ■ v M = e M u M 
-B ■ u M - A ■ v M = e M v M . (4) 

Here A and B are real L x L matrices whose non-zero 
elements are given by A^, = —Thi, Aj^+i = Aj+i^ = 
-Ji/2, = -B i+ li = -Ji/2. ' (For the PBC 

spin-chain case, we have the additional matrix elements 
A LA = A hL = (J L /2)(-l) N *, and B LA = -B ljL = 
(J l /2)(-1) Nf ). While in the ordered case the solution 
of Eqs.((3]) can be reduced, by switching to momentum- 
space, to independent 2x2 problems, in the general dis- 
ordered case one has to diagonalize the 2L x 2L problem 
numericall y 37 ' 38 . 

The spectrum of Eqs. 0$ turns out to be given by ±e M , 
with e M > 0, and in terms of the new fermion operators, 
H becomes: 

L L ^ 

^i— 1 M— 1 

The ground state of H is the Bogoliubov vacuum state 
|\&o) annihilated by all 7^ for /j, = 1---L, 7 jU |^'o) = 0, 
with an energy E = - X^=i 

2. Dynamics 

The Schrodinger dynamics associated to a time- 
dependent H(t) can be solved by a time-dependent Bo- 
goliubov theory 3 ^. The basic fact that makes the solu- 
tion possible even in the time-dependent case is that the 
Heisenberg's equations of motion for the operators Ci t H{t) 
are linear, because the Hamiltonian is quadratic: 

d L 
ih-c^ H {t) = 2^ [A^{t)c 3 , H {t) + B^(i)c] H {t) . 

i=i 

(6) 

Here the matrices A and B have the same form given 
previously, except that now the time-dependence of T(t) 
is explicitly accounted for. If we denote by 7M-i R 

the Bo- 
goliubov operators that diagonalize H(t- m ) at the initial 
time, and ujf , v™ the corresponding initial eigenvectors, 
it is simple to verify that the Ansatz 

L 

CiM t ) = ^( U ip( t h^ + <»( t hlin) > (7) 

11=1 
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does indeed solve the Heisenberg equations ((6|) , provided 
the time-dependent coefficients Ui a {t) and v^it), satisfy 
the following system of first-order differential equations 

d 2 L 

V*® = ~\^ ^'^>^ + -(8) 

with initial condition u iu (t in ) = uj", v ifi (t in ) = v™. 
Eqs. ([8]) are the natural time-dependent generalizations 
of the static Bogoliubov-de Gennes Eqs. ||U), and, once 
again, they have to be solved numerically in the general 
disordered case. 



III. RESIDUAL ENERGY AND KINK DENSITY 



7/j,inj 7u in: u i,fi{^)i an d ^i,/i(0), and uses the fact that the 
lu.in annihilates by construction the initial state |\l/(ti n )). 

By applying this procedure to the calculation of Ea n 
we get: 

E Rn = £ (^(0) H0)v\0) - u(0)«t( )] y + 

id 

BiM[v(Q)u1(0)-u(0)vt(pj\.?) ,(11) 

where u(0) and v(0) are L x L matrices with elements 
Ui tfl (0) and Vi^(0). Similarly, the density of defects pk 
can also be expressed as: 

P* = ^L%{ 1 - (K0)-<0)] [nt(0)+«t(0)]). . 

(12) 



How effectively the Schrodinger dynamics drives the 
system from the initial disordered quantum ground state 
l^in) towards the classical ground state \^{) = Yli IT)i 
(or the fully reversed one = 14)0? 

A way of quantifying the degree of adiabaticity of the 
evolution is given by the residual energy, defined as 

-E^res = ^fin — ^cl j (9) 

where E c \ = — J2i Ji i s the classical ground state energy 
of H(t Rn ) = H ch and £ fin = (*fi„|F c i|*fin) is the aver- 
age classical energy of the final time-evolved state l^fin). 
Obviously, E^ ni and hence -E ros , depends on the param- 
eters specifying the evolution: the smaller and closer to 
E c \ the "slower" the evolution. 

An alternative way of quantifying the degree of adia- 
baticity of the evolution, is given in terms of the density 
of kinks pk in the final state, defined by 

pu = ^X>(o)l5 (i-^f+i) |*(o)> (io) 

i 

(for a PBC chain the sum goes up to L, instead of L— 1). 

When no disorder is present the two quantities coin- 
cide, apart from trivial constants. In the disordered case, 
however, this is not the case. A defect will form with 
higher probability at a link where the corresponding ex- 
change coupling Jj is small. Therefore the residual en- 
ergy is not simply given by the kink density times the 
exchange coupling. 

The calculation of quantities like E^ n or pk is 
straightforward. Quite generally, given an opera- 
tor 0[cj,cj] expressed in terms of the c,'s and cj's, 
its expectation value over the final state I* (tan = 
0)) can be expressed, switching from the Schrodinger 
to the Heisenberg picture, as (*(0)|O[ci, cj]|*(0)) = 
(*(i in )|O[c ii// (0),c l t H (0)]|*(ii n )). Next, one uses the ex- 
pressions ((Tj) for the c^h (0)'s and c\ H (0) in terms of 



IV. RESULTS 

Our results for the dynamics are obtained by inte- 
grating numerically the time-dependent Bogoliubov-de 
Gennes equations ([8]). As initial point of the evolution 
it is enough to consider t{ n — — 5t, taking u™ and v™ 
from the diagonalization of H{t- m ) according to Eq. (Qj: 
we checked that our results do not depend on the precise 
value of ti n , as long as it is not too small. We consid- 
ered systems up to L — 512 and annealing times up to 
t = 1000. Ensemble averages are calculated over a suit- 
ably large number of disorder realizations (of the order 
of 1000). The analysis of the instantaneous spectrum 
and its statistics has been obtained by solving the static 
Bogoliubov-de Gennes eigenvalue equations ([!]) for sys- 
tems up to L — 512. 



A. Landau-Zener transitions and disorder 

In order to get an initial understanding on the mecha- 
nisms that lead to breaking of adiabaticity in the present 
system, it is instructive to consider in more detail the 
time-evolution of a single realization of the disorder. To 
be specific, Fig. fl} shows the time-evolution of the resid- 
ual energy E les (t) = (#(t)|JT(t)|¥(t)) - E gs (T(t)), where 
£7 gs (r(t)) is the instantaneous ground state energy cor- 
responding to r(t), for a single L = 64 sample and for 
values of r up to 5000. We also plot the instantaneous 
spectral gaps of the problem (thick solid lines) obtained 
by diagonalizing the Hamiltonian for any given value of 
the parameter T. As mentioned previously, the dynamics 
conserves the fermion parity, so that only excitations in 
the same fermion parity sector are accessible. If we order 
the single-particle eigenvalues as ei < ti < ■ ■ ■ < ej,, then 
the lowest excited state accessible to the dynamics (i.e., 
conserving the fermionic parity) is associated with an 
excitation energy Ai = 2(ei + 62), rather than A = 2e%. 
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FIG. 1: (Color online) Residual energy E les (t) versus t for a 
given instance with L — 64 of the random Ising model with 
transverse field, at different values of r. The solid lines are 
the lowest-lying instantaneous spectral gaps A n as a function 
of r. 



The next excited state is A2 = 2(ei + 63), and so on. 
These are the instantaneous gaps shown in Fig. (JTj). 

An important feature which emerges from this exam- 
ple is that one cannot in general locate a single specific 
value of r where the minimum and most important gap 
is present. Certainly, typically the first occurrence of a 
small gap during the annealing trajectory is close to the 
critical point, T c = 1. Usually, this critical-point gap is 
also the smallest one that the systems encounters dur- 
ing its evolution. However, it can happen, as Fig. ([T]) 
shows, that the system safely goes through the critical- 
point small gap (see t = 5000 results) but then looses 
adiabaticity due to a comparable gap encountered later 
on (here at T ~ 0.5). Once adiabaticity is lost, the sys- 
tem will generally miss to follow the first excited state 
either, getting more and more excited as time goes by. 

It seems clear that the analysis of the adiabatic dy- 
namics of a disordered Ising chain requires a knowledge 
of the statistics of these low-lying gaps in the spectrum 
(in the pertinent parity sector). We concentrate our at- 
tention on the region close to the critical point, where 
the smallest gaps are found, for large L. 

We start asking how these smallest gaps are dis- 
tributed, for different realizations of the disorder. Let 
us denote by P(Ai,L) the distribution of gaps Ai = 
2(ei + 62) (the lowest one relevant for the dynamics) 
for a chain of length L, assumed to be normalized: 
/ °°dAi P(Ai,L) = 1. For the smallest gap A = 2ei, 
Young and Rieger— have shown that the correct scal- 
ing variable which makes the critical point distribution 
universal, for different L, is — log(A)/vL- By using a 
scaling variable of the same form, g = — log(Ai)/\/Z, 
we see that the gaps Ai are also distributed in the same 
universal way, see Fig. This implies that at the crit- 
ical point, P*(g) = VTe-fvTp^-gvT. £) ^ for large 
L, universal and normalized. As a consequence, gaps at 
the critical point have an extremely wide distribution, 
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FIG. 2: (Color online) Distribution of Ai = 2(ei + e 2 ), 
the smallest gap relevant for the dynamics, at the critical 
point r c = 1 for different systems sizes, showing the col- 
lapse of the distributions P(Ai, L) when the scaling variable 
g — — log (Ai)/VX is used. The resulting distribution is the 
P*(g) discussed in the text. 



for large L, with typical gaps which are exponentially 
small^^^ in the system size: [Ai] typ oc e _Cv ^. 

B. Density of kinks 

Given the wide distribution of the instantaneous gaps, 
it is important to understand how this reflects itself in the 
distribution of various observables. We first consider the 
behavior of the density of defects pk defined in Eq. (fT0|) . 
The results for the probability distribution function of 
Pk, P{Pk), ar e presented in Fig. (J3j) for r = 10 and r = 
1000. The distribution P(pk), for given r, is found to be 
approximately log-normal: 



PiPk) 



1 



_L e -(tePk-lnpk)' d /2<rl 



with a standard deviation ol decreasing as 1/yL. The 
data collapse of the results for different L, in terms of 
the variable (hapk — hipk)/crL, shown in the inset, qual- 
ifies the accuracy of this statement. This -\/i-reduction 
of the width of the log-normal distribution P{pk) with 
increasing L is at variance with the result obtained for 
the distribution of the gaps at the critical point, whose 
width increases as \[L; here, on the contrary, the correct 
scaling variable appears to be (In pk — ]&pk)"\/L, rather 
than (lnpfc — In pk)j\fL. This width reduction, for in- 
creasing L, implies that the average density of defects 
[pk]av approaches the typical value [pfe]typ = e' lnpfc l av for 
large enough L, since [pk\&v = e lnpfc+<T ^/ 2 implies that: 



[p 



fcjav 



[p. 



fejtyp 



[Pk. 



typ 



= e oi/2 _ 1 * 

L 



(13) 



This fact is shown explicitly in Fig.lfJJ (top), where we 
see that large deviations between [/Ofc]t yp = eP npfc l av and 
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FIG. 3: (Color online) Probability distribution for the loga- 
rithm of the density of defects x = — In pt , for two different 
annealing rates r. The distribution function is universal and 
log-normal with a variance ol which scales as l/y/T. In the 
insets we show the data collapse of all the curves when plot- 
ted as a function of the reduced variable (x — x)/<tl, where 
x = — In pk . 



[pk]av are seen only for L < 64. For large systems, 
L > 128, the two quantities are essentially coincident, 
for all values of r. Despite the universal behavior of the 
distribution P(pk) at all annealing rates, the behavior 
of [/3fe] av (T) changes drastically between short and long 
T >g33 pjg (gj (bottom) focuses on the average kink den- 
sity [p/c]av for various L, as a function of r. The initial 
small-r behavior of [pk]&v(i~), indicated by the dashed 
line in Fig. ([4]), seems a power-law, [pk]a,v(T~) ~ t~ - 5 , 
i.e., exactly what one finds for the ordered Ising chain 8 , 
where the result is interpreted in terms of the the Kibble- 
Zurek mechanism. A possible explanation resides in the 
fact that our model presents a Griffiths phase extending 
for all r > r<4£. This phase is characterized by a gap 
A ~ L~ z , where the dynamical exponent z(T) is a contin- 
uous function of the parameter T, diverging, z — > oo, for 
r — > r c , while saturating to a constant for large T. The 
second gap, which is relevant for our dynamical problem, 
shows a similar behavior^ Ai ~ L~ z , with a dynam- 
ical exponent z'(T) = z(T)/2. For fast annealing rates, 
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FIG. 4: (Color online) Top: Comparison between average 
[/9fc]av and typical [pfc]t yP = e l ln Pfc]=v k m k density for differ- 
ent system sizes on varying the annealing rate r. The same 
symbol is used for both cases. The typical value (dashed line) 
lies always below the average value (continuous line) , but the 
difference between the two is negligible for L > 128. Bottom: 
Average kink density [p^Jav as a function of the annealing rate 
r for chains of different lengths L = 16, 32, 64, 128, 256, 512. 
The data for [pfe] a v are the same appearing in the top part of 
the figure. The dashed line is a power-law describing the 
small-r behavior, [pk\a.v(r) ~ r _0 j . The solid thick line 
through the [pfc] a v data is a fit with a function A/ log 2 (jr), 
described in the text. The averages are calculated over 1000 
different realizations of disorder. 



the system loses adiabaticity before reaching the criti- 
cal point, well inside the T > r c Griffiths phase. As 
in the ordered case, the gaps exhibited by such a phase 
would induce a defect density decreasing as a power-law 
of the annealing time t, with the crucial difference that 
the power-law exponent is not constant here, due to the 
T-dependence of z' . One should expect, presumably, a 
gradual crossover with a power-law exponent which be- 
comes smaller and smaller, connecting in a gentle way 
with the large r behavior of [pfc] a m which shows marked 
deviations from a power-law behavior. Dziarmaga, based 
on scaling arguments 33 showed that at large r the den- 
sity of kinks should decrease as the inverse square of the 
logarithm of r. Our data for the largest systems agree 
very well with this prediction, as the best-fit (solid line 
in Fig. ([?])) shows. 

A bound to [pk]av{T) can also be constructed by a 
Landau-Zener argument — complemented by a knowl- 



edge of the distribution of the first gap P(Ai,L) — , in 
a similar fashion to that presented by Zurek et al& for 
the ordered Ising case. The derivation starts by con- 
sidering the probability P cx (r, L) of loosing adiabaticity 
for a system of size L, when turning off Y with an an- 
nealing rate r" 1 . Evidently, P cx (t,L) > P e c ^- point (r, L), 
where we have denoted by P e "' pomt (r, L) the probability 
of getting excited by Landau-Zener events at the criti- 
cal point (indeed, we have seen that there is a chance 
of getting excited also by gaps well below the critical 
point). P C x' point (r, L), in turn, can be constructed by 
knowing the distribution of the gaps Ai at the critical 
point, and the simple two- level Landau-Zener formula 



3 LZ 



(a being the slope of the two ap- 



proaching eigenvalues). Lumping all constants together, 
7 = 7r/(4/to), we write P^ = e nAl and assume that 
the distribution of 7 cx a -1 is not important in our esti- 
mate, while that of Aj is, so that: 



P«-P° int (r,L) = ^ d^P{^ u L)e 
Jo 

dg P*(g) e~^ 2VTg , (14) 



-1\f~Lg 



where the second equality follows from switching to the 
scaling variable g = — log (Ai)/\/Z. Obviously, for r = 
we correctly have P c c x rpoint (r = 0,L) = dg P*(g) = 
1, from the normalization condition. When r is finite, 

the LZ factor e~ 7Te La provides a lower cut-off in the 
integral at a characteristic g c = log (jr) / (2\/~L) , and this 
cut-off is sharper and sharper as L increases: one can 

verify that, for large L, e~ 1Te 2 ^ Za w 9(g — g c ). 
consequence, for large enough L we can rewrite: 



As a 



i£- potat (T,£)«nCff C 



dgP*{g), (15) 



i.e., P^' point (T, L) turns out to be a universal function 
of the scaling variable g c = log (jT)/(2y/L), for L large. 
This universal function H(g c ) is shown in Fig.©, where 
we see that data for L > 512 collapse into a single curve. 
The density of kinks for large r, and large enough L, 
can be obtained by evaluating the typical length L e (r) 
of a defect-free region upon annealing, e being a small 
quantity of our choice, denoting the probability of getting 
excited. Since P™ pomt (r, L) w YV{g c ) is a lower bound for 
P ex (T, L), we have that 



(16) 



where LI -1 denotes the inverse function of II. If we 
now identify the inverse of the defect-free region length, 
I/7 1 (r), with the density of kinks /3fc(r), we get the fol- 
lowing lower bound for the latter: 



Pk{r) 



1 



L e (r) ~ log 2 ( 7 r) 



> 



(17) 
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FIG. 5: (Color online) Approach to the universal function 
n(<? c ) for increasing chain lengths L, see text. All data from 
L > 512 collapse well into a single curve. Inset: P c C x' p ° lnt (r, L) 
obtained from the integral in Eq. (|14|l versus r for different 
values of L. 



On the basis of this argument, we conclude that the den- 
sity of kinks cannot decrease faster than 1 / log 2 (77-) for 
large r, which agrees with the argument discussed by 
Dziarmaga^. 



C. Residual energy 

In the ordered case the residual energy per spin is 
simply proportional to the kink-density, E res / L — 2Jpk, 
while here, evidently, kinks sitting at small J^'s are fa- 
vored, on average, by the adiabatic evolution process. It 
is therefore of importance to analyze the scaling of the 
residual energy that, as we will show, differs quantita- 
tively from that of the kink density. Since kinks will be 
formed on the weak links, one expects on general grounds 
that the residual energy would decay faster than the kink- 
density for large r's. 

As in the case of the kink density, we first analyze 
the probability distribution for the residual energy per 
site, which we present in Fig.©. Once again the resid- 
ual energies are approximately log-normal distributed 
and can be reduced to a universal form (see the in- 
sets) when properly rescaled, i.e., in terms of the variable 
(In (E Ies /L) - In (E les /L))VL. 

The average residual energy per site [E res /L] av as a 
function of the annealing time r shows a crossover from 
a power-law decay, approximately t -1 for fast quenches, 
to a much slower decay (see below) for slow evolutions. It 
is interesting to note that although for fast quenches the 
disorder is considered to play a minor role, nevertheless 
the exponent of the decay of the residual energy differs 
from that of the kink density. The analysis of the regimes 
of large r's is more delicate. The LZ argument given 
above tells us nothing about the behavior of the residual 
energy for large r. We then proceed as follows. Assuming 
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FIG. 6: (Color online) Probability distribution for the resid- 
ual energy per site at two different annealing rates r _1 . The 
distribution function is universal and log-normal with a vari- 
ance which scales as \/y/~L. In the insets we show the data 
collapse. 
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FIG. 7: (Color online) Top: Average residual energy 
per site [E lea / 'L] av as functions of the annealing rate r for 
chains of different lengths L = 16, 32, 64, 128, 256, 512. The 
dashed line is the power-law describing the small-r behavior, 
[E TeB / 'L]a V (r) ~ r _1 . Averages are calculated over 1000 re- 
alizations of disorder. Bottom: The ratio of the density of 
kinks and the residual energy versus r, used to extract the 
power of the log- dependence of E rea . 



for the residual energy a logarithmic behavior similar to 
that found for pk 



1 



log c (7T) ' 



(18) 



we can determine £ from the data of Fig. (©(Top) by 
plotting the ratio of [pfe]"v an d [-E'res/^Jav versus t for 
several values of a, as done in Fig. (©(Bottom). If 
[pk]ov ~ log 2 il T )i then the value of a which makes 
this ratio constant is: 



[E ies /L] 



oc log^ 2a (77-) ~ const. 



C/2 .(19) 



Numerically, see Fig. 0, we find a w 1.7 ± 0.1, which 
implies ( ~ 3.4 ± 0.2. 



V. DISCUSSION AND CONCLUSIONS 

In this paper we have studied the adiabatic quantum 
dynamics of a one-dimensional disordered Ising model 



across its quantum critical point. Our main results can 
be summarized in the dependence of the average kink 
density (see however—) and residual energies as a func- 
tion of the annealing rate 



[E T , 



[pk 
[pk 



r -0.5 

(logr) 



(logr) 



fast quenches 
large r , 

fast quenches 
large r, with £ <• 



3.4 



Although the dynamics is dominated by a very wide dis- 
tribution of gaps at the critical point, P*(— In (Ai)/y/~L) 
(see Fig. ((5J)), we find that the distribution for both these 
quantities are log-normal but with a variance that de- 
crease, like 1/VZ, for increasing chain length L: typical 
and average values, therefore, coincide for large L. The 
wide distribution of gaps, on the other hand, with its 
characteristic ln(Ai)/vL scaling, is responsible, within 
a Landau-Zener theory, for the extremely slow decay of 
the average density of kinks, [pfc] av ~ l/(lnr) 2 . This dis- 
cussion applies only for reasonably large sizes L. If L is 



9 



small, the minimum gap Ai of a given instance can be 
sufficiently large that the adiabatic regime, predicted to 
occur beyond a characteristic r c oc A^~ 2 , is actually seen: 
a fast decay of and E ICS /L is expected 41 for t > r c , in 
such a case. 

It is interesting to compare these results with those of 
a classical thermal annealing, where, according to Huse 
and Fisher—, the residual energy also shows a logarithmic 
behavior 

^(r)/L~(logr)- CoA Cca<2, 

but with an exponent Cca which is bound by C,ca < 2. 

If we look at this problem from the perspective of op- 
timization algorithms, it seems that quantum annealing 
(QA) gives a quantitative improvement over classical an- 
nealing for the present system, as is indeed found in other 
cases^Mi^i 4 ^! 4 ^, but not always (Boolean Satis- 
fiability problems seem to be a test case where QA per- 
forms worse than classical annealing, see Ref. 50). 

In this respect, however, several important issues re- 
main to be clarified. First of all, AQC-QA has a large 
freedom in its construction: the choice of the possible 
source of quantum fluctuations 5 - 1 - — generally speaking, 
one can take H(t) = Ha n + J2\F\(t)H\ — , and the 
time-dependence of the various T\ (t) , which need not be 
linear in tim o 12 ' 13 . Regarding the time dependence of 
the couplings, we simply note that an optimal choice of 
the "speed" T(t) with which the critical point is crossed 
can provide an improvement in the exponents 1 - 3 -, but def- 
initely not change a logarithm into a power-law. Regard- 
ing the possibility of adding extra kinetic terms to H(t), 
it is clear that terms like — T xy (t) JY J%o\ °f+i (XY- 
anisotropy) or similar short range interactions will not 



change the universality class of the infinite randomness 
quantum critical point of the present model 3 - 4 -. Hence, a 
logarithmically-slow AQC-QA is expected also in more 
general circumstances, for the present one-dimensional 
model. We expect this to be a genuine consequence of 
the randomness present in the problem at hand, which 
makes the adiabatic quantum dynamics intrinsically slow 
and ineffective in reaching the simple classical ferromag- 
netic ground state o 52 ' 53 . 

This is perhaps to be expected in view of the results 
of VidalS 4 -, who showed that problems where the entan- 
glement entropy of a block is bound, can be computed 
classically with a comparable efficiency. Generically, in 
disordered one-dimensional system the entanglement en- 
tropy grows at most logarithmically with the system size 
at a critical poin t 55 ' 56 ! 57 , at this is not enough to sub- 
stantially change the relative efficiency of quantum ver- 
sus classical algorithms. 

Therefore, the route to investigate seems to be follow- 
ing: search for models in more then one-dimension, where 
the entropy of entanglement grows stronger, which, at 
the same time, have "gentle" enough critical point gap 
distributions. 
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